#Load data
GetData<-function(size=5000,boot=TRUE,directory="",y0=1986,y1=2015){
  if(directory==""){
    Datawd<-getwd()
  }else{
    Datawd<-directory
  }
  yy0<-substr(y0,(nchar(y0)+1)-2,nchar(y0))
  yy1<-substr(y1,(nchar(y1)+1)-2,nchar(y1))
  #Load data for years y0 and y1
  load(file.path(Datawd,paste0("Data",yy0,"-",yy1,".Rda")))
  
  vars<-c("lrwage3","union","pubsect","manuf","nonwhite",
          "female","partte","married","smsa",
          "exper" , "exper2" ,
          "exp_cl" , 
          paste0("expe",1:length(unique(Data$exp_cl))), 
          "years_educ" , "educ_cl" , 
          paste0("educ",1:length(unique(Data$educ_cl))), 
          "edex" ,
          "region",
          paste0("reg",1:length(unique(Data$region))),
          "ee_cl",
          paste0("ee",1:length(unique(Data$ee_cl))),
          "ind21",
          paste0("indu",1:length(unique(Data$ind21))),
          "state",
          paste0("state",1:length(unique(Data$state))),
          "orgwgt","year")
  
  Z0<-Data[Data$year==y0,vars]
  Z1<-Data[Data$year==y1,vars]
  if (boot==TRUE){
    w="orgwgt"
    probab0=eval(parse(text=paste("Z0$",w,"/sum(Z0$",w,")",sep="")))
    set.seed(1970)
    i0 <- sample(x=1:nrow(Z0),size=size, replace = FALSE,prob=probab0)
    Z0<-Z0[i0,vars]
    
    probab1=eval(parse(text=paste("Z1$",w,"/sum(Z1$",w,")",sep="")))
    set.seed(1983)
    i1 <- sample(x=1:nrow(Z1),size=size, replace = FALSE,prob=probab1)
    Z1<-Z1[i1,vars]
  }
  
  rm(Data)
  
  return(list(Z0=Z0,Z1=Z1))
}

#Function to get coefficients from the QR
ParRq<-function(t,f=formula,d=Z0,w="orgwgt",mth="pfn"){
  require(quantreg)
  if (w!=""){
    h<-eval(parse(text=paste("rq(f, data=d,tau=t,weights=",w,",method='",mth,"')$coef",sep="")))
  } else{
    h<-rq(f, data=d,tau=t,method=mth)$coef
  }
  return(h)
}

#Function to get the coefficients from the QR in parallel
ParCoefRq<-function(J,formula,Data,w="orgwgt",mth="pfn",nodes=13,inorder=T,ParRq){
  require(quantreg)
  require(doParallel)
  taus <- 1:J/(J+1)
  if(nodes>1){
    cl <- makeCluster(nodes, type = "SOCK")
    registerDoParallel(cl)
    CoefRq<-foreach(tt=taus,.combine='cbind',.inorder=inorder)%dopar%{
      tryCatch(
        {
          ParRq(tt,f=formula,d=Data,w=w,mth=mth) 
        },
        error=function(cond){
          set.seed(13)
          if (w!=""){
            probab=eval(parse(text=paste("Data$",w,"/sum(Data$",w,")",sep="")))
            iboot <- sample(1:nrow(Data), replace = TRUE,prob=probab)
          }else{
            iboot <- sample(1:nrow(Data), replace = TRUE)
          }
          DataB <- Data[iboot,]  
          coe<-ParRq(t,f=formula,d=DataB,w=w,mth=mth)
          return(coe)
        }
      )
    }
    stopCluster(cl)
  }else{
    CoefRq<-NULL
    for (tt in taus) {
      CoefRq<-cbind(CoefRq,ParRq(tt,f=formula,d=Data,w=w,mth=mth))
    }
  }
  return(CoefRq)
}

Impact<-function(n=6,maxIter=10,seed=T,w="",Q=c(0,0.25,0.5,0.75,1),
                 mth="pfn",nodes=1,J,s,formula,h,Data,
                 ParRq,Legendre,polynomial.values){
  require(quantreg)
  require(orthopolynom)
  
  p<-length(h)
  q<-length(Q)
  
  if(seed==T){set.seed(s)}
  if (w!=""){
    probab<-eval(parse(text=paste("Data$",w,"/sum(Data$",w,")",sep="")))
    iboot <- sample(1:nrow(Data), replace = TRUE,prob=probab)
  }else{
    iboot <- sample(1:nrow(Data), replace = TRUE)
  }
  DataB <- Data[iboot,]  
  taus <- 1:J/(J+1)
  ##Sometimes the DataB produces error in the rq
  CoefRq<-ParCoefRq(J,formula=formula,Data=DataB,w=w,mth=mth,nodes=nodes,ParRq=ParRq)
  
  IIn<-matrix(rep(0,q*p),ncol=p,nrow=q)
  for (pp in 1:p){
    eval(parse(text=paste0('ft<-CoefRq[',pp,',]')))
    for (qq in 0:(q-2)){
      if (qq==0){
        q_fin=Q[qq+2]
        q_ini=0
      }else if (qq==(q-2)){
        q_fin=1
        q_ini=Q[qq+1]
      }else{
        q_fin=Q[qq+2]
        q_ini=Q[qq+1]
      }
      IIn[qq+1,pp]<-II(ft,taus,q_ini,q_fin,n)
    }
    IIn[q,pp]<-II(ft,taus,0,1,n)
  }
  
  colnames(IIn)<-h
  return(2*IIn)
}

polynomial.values <- function(polynomials, x, matrix=FALSE){
  ## Changed copy of polynomial.vales from orthopolynom in order
  ## to add matrix argument
  require(polynom)
  n <- length(polynomials)
  if(!matrix) {
    values <- vector(mode="list", length=n)
  } else {
    values <- matrix(ncol=n, nrow=length(x))
  }
  j <- 1
  while(j <= n) {
    if(!matrix) {
      values[[j]] <- predict(polynomials[[j]], x)
    } else {
      values[, j] <- predict(polynomials[[j]], x)
    }
    j <- j + 1
  }
  values
}

Legendre <- function(x, n, normalized=TRUE, intercept=FALSE){
  ## Create a design matrix for Legendre polynomials
  ## x - numeric
  ## n - see orthopolynom
  ## normalized - logical, see orthopolynom
  ## intercept - logical, add intercept
  tmp <- legendre.polynomials(n=n, normalized=normalized)
  if(!intercept) tmp <- tmp[2:length(tmp)]
  polynomial.values(polynomials=tmp, x=x, matrix=TRUE)
}

II_L<- function(x, n, normalized=TRUE, intercept=FALSE){
  tmp <- legendre.polynomials(n=n, normalized=normalized)
  if(!intercept) tmp <- tmp[2:length(tmp)]
  tmp <- polynomial.integrals(tmp)
  tmp <- polynomial.integrals(tmp)
  polynomial.values(polynomials=tmp, x=x, matrix=TRUE)
}

II<-function(ft,taus,a,b,n){
  Mod<-lm(ft ~ Legendre(taus,n,normalized=TRUE,intercept=FALSE))
  k<-Mod$coefficients[names(Mod$coefficients)=="(Intercept)"]
  v<-(II_L(b,n)-II_L(a,n))%*% Mod$coefficients[names(Mod$coefficients)!="(Intercept)"]+0.5*(k)*(b^2 - a^2)
  return(v)
}

#variance by column
colVars <- function(x, na.rm=FALSE, dims=1, unbiased=TRUE, SumSquares=FALSE,
                    twopass=FALSE) {
  if (SumSquares) return(colSums(x^2, na.rm, dims))
  N <- colSums(!is.na(x), FALSE, dims)
  Nm1 <- if (unbiased) N-1 else N
  if (twopass) {x <- if (dims==length(dim(x))) x - mean(x, na.rm=na.rm) else
    sweep(x, (dims+1):length(dim(x)), colMeans(x,na.rm,dims))}
  (colSums(x^2, na.rm, dims) - colSums(x, na.rm, dims)^2/N) / Nm1
}
w_mean<-function(x,prob){
  return(sum(x*prob))
}

change<-function(W0,W1,Q=c(0,0.25,0.5,0.75,1),seed=69,m=4500,Boot=TRUE){
  if(Boot==TRUE){
    w0<-BootW(W0,m=m,seed=seed)
    w1<-BootW(W1,m=m,seed=seed)
  }else{
    w0<-W0
    w1<-W1
  }
  
  M0<-Gini(w0,n=6,Q=Q)
  M1<-Gini(w1,n=6,Q=Q)
  return(M1-M0)
}

changeG<-function(Z0,II0,Z1,II1,h,coe,Q=c(0,0.25,0.5,0.75,1),seed=69,m=4500,Boot=TRUE){
  firstup <- function(x) {
    substr(x, 1, 1) <- toupper(substr(x, 1, 1))
    x
  }
  if(Boot==TRUE){
    z0<-BootZ(Z0,m=m,seed=seed)
    z1<-BootZ(Z1,m=m,seed=seed)
  }else{
    z0<-Z0
    z1<-Z1
  }
  X0_bar<-X_bar(z0,h=h)
  X1_bar<-X_bar(z1,h=h)
  
  Gini0_hat<-Gini_hat(X0_bar,II0)
  Gini1_hat<-Gini_hat(X1_bar,II1)
  Gini10_hat<-Gini_hat(X0_bar,II1)
  Coef<-Gini10_hat-Gini0_hat
  Cova<-Gini1_hat-Gini10_hat
  
  ### Impacts
  for (cc in coe) {
    eval(parse(text=paste0('X1_bar_',cc,'0<-Ind_cov(X0_bar,X1_bar,"',cc,'")')))
    eval(parse(text=paste0('100*II1["',cc,'",]*(X0_bar["',cc,'"]-X1_bar["',cc,'"])')))
    eval(parse(text=paste0('Gini10_hat_',cc,'<-Gini_hat(X1_bar_',cc,'0,II1)')))
    eval(parse(text=paste0(firstup(cc),'<-Gini1_hat-Gini10_hat_',cc)))
    
    eval(parse(text=paste0('II0_',cc,'1<-Ind_coef(mean0,mean1,"',cc,'")')))
    eval(parse(text=paste0('Gini0_',cc,'1<-Gini_hat(X0_bar,II0_',cc,'1)')))
    eval(parse(text=paste0(firstup(cc),'C<-Gini0_',cc,'1-Gini0_hat')))
  }
  
  eval(parse(text=paste0('O<-matrix(c(Cova,Coef,',paste(firstup(coe),collapse=','),',',
                         paste(paste0(firstup(coe),'C'),collapse=','),'),nrow=length(Cova),ncol=',2+2*length(coe),')')))
  
  eval(parse(text=paste0('colnames(O)<-c("Cova","Coef",',paste(paste0('"',firstup(coe),'"'),collapse=','),',',
                         paste0(paste0('"',paste0(firstup(coe),'C'),'"'),collapse=','),')')))
  return(O)
}

#Function get bootstrap confidence interval in paralel
ParBootCI<-function(R=1000,nodes=15,
                    BootW,Gini,Lorenz,BootZ,X_bar,Gini_hat,Ind_cov,Ind_coef,
                    mean0,mean1,
                    fun=change,...){
  require(doParallel)
  cl <- makeCluster(nodes, type = "SOCK")
  registerDoParallel(cl)
  
  Boot<-foreach(r=1:R,.combine='cbind',.inorder=FALSE)%dopar%{
    BootW<-BootW
    Gini<-Gini
    Lorenz<-Lorenz
    BootZ<-BootZ
    X_bar<-X_bar
    Gini_hat<-Gini_hat
    Ind_cov<-Ind_cov
    Ind_coef<-Ind_coef
    mean0<-mean0
    mean1<-mean1
    fun(...,seed=r)
  }
  stopCluster(cl)
  return(Boot)
}

BootW<-function(W,m,seed=69){
  set.seed(seed)
  ii<-sample(x=1:length(W),size=m, replace = TRUE)
  Bw<-W[ii]
  return(Bw)
}

Gini<-function(Z,n=6,Q=c(0,0.25,0.5,0.75,1)){
  library("orthopolynom")
  q=length(Q)
  LC<-Lorenz(Z,fig=FALSE)
  leg_fun<-polynomial.functions(legendre.polynomials(n=n, normalized=TRUE))
  x <- -1+2*LC$tau
  p.list<-as.list(rep(0,length(leg_fun)))
  for (i in 1:length(leg_fun)){
    p.list[[i]]=as.function(alist(x = , leg_fun =  , leg_fun(2*x-1)))
  }
  
  X=''
  for (j in 0:n){
    v=paste("p.list[[",j+1,"]](LC$tau,leg_fun[[",j+1,"]])",sep="")
    X=paste(X,v,sep=",")
  }
  X=substr(X,2,nchar(X))
  X=paste("cbind(",X,")")
  X=eval(parse(text=X))
  
  p.list<-legendre.polynomials(n=n, normalized=TRUE)
  p.int <- polynomial.integrals( p.list )
  
  #First Integral of the legendre polynomials
  I_f<-polynomial.functions(polynomial.integrals(p.list))
  I_p<-as.list(rep(0,length(I_f)))
  for (i in 1:length(leg_fun)){
    I_p[[i]]=as.function(alist(x = , I_f =  , 0.5*I_f(x)))
  }
  I<-matrix(rep(0,q*(n+1)),ncol=n+1,nrow=q)
  
  for (qq in 0:(q-2)){
    if (qq==0){
      q_fin=Q[qq+2]
      q_ini=0
    }else if (qq==(q-2)){
      q_fin=1
      q_ini=Q[qq+1]
    }else{
      q_fin=Q[qq+2]
      q_ini=Q[qq+1]
    }
    TI=''
    for (j in 0:n){
      ii=paste("I_p[[",j+1,"]](",2*q_fin-1,",I_f[[",j+1,"]])-I_p[[",j+1,"]](",2*q_ini-1,",I_f[[",j+1,"]])",sep="")
      TI=paste(TI,ii,sep=",")
    }
    TI=substr(TI,2,nchar(TI))
    TI=paste("cbind(",TI,")")
    I[qq+1,]=eval(parse(text=TI))
  }
  q_fin=1
  q_ini=0
  TI=''
  for (j in 0:n){
    ii=paste("I_p[[",j+1,"]](",2*q_fin-1,",I_f[[",j+1,"]])-I_p[[",j+1,"]](",2*q_ini-1,",I_f[[",j+1,"]])",sep="")
    TI=paste(TI,ii,sep=",")
  }
  TI=substr(TI,2,nchar(TI))
  TI=paste("cbind(",TI,")")
  I[q,]=eval(parse(text=TI))
  
  Mod<-lm(LC$L ~ X-1)
  
  AreaUndL<-apply(apply(I,1,"*",Mod$coef),2,sum)
  Gini<-matrix(rep(0,q),ncol=1,nrow=q)
  for (qq in 0:(q-2)){
    Gini[qq+1]<-(Q[qq+2]-Q[qq+1])*Q[qq+1]+0.5*(Q[qq+2]-Q[qq+1])^2-AreaUndL[qq+1]
    Gini[qq+1]<-2*Gini[qq+1]
  }
  Gini[q]<-0.5-AreaUndL[q]
  Gini[q]<-2*Gini[q]
  return(Gini)
}

Lorenz<-function(Z,fig=TRUE){
  L<-cumsum(sort(Z))/sum(Z)
  tau<-1:length(L)/length(L)
  if (fig==TRUE){
    plot(tau,L,type="l",lwd = 2,xlim=c(0,1),ylim=c(0,1),xlab="Cumulative share of people from lowest to highest", ylab="Cumulative share",xaxs="i", yaxs="i")
    lines(tau,tau,lwd=2)
    polygon(c(tau,rev(tau)),c(tau,rev(L)),col="mediumseagreen")
  }
  return(list(tau=tau,L=L))
}

BootZ<-function(Z,m,seed=69){
  set.seed(seed)
  ii<-sample(x=1:NROW(Z),size=m, replace = TRUE)
  Bz<-Z[ii,]
  return(Bz)
}

X_bar<-function(X,h,w="orgwgt"){
  require(matrixStats)
  prob<-t(as.matrix(X[w]/sum(X[w])))
  X<-as.matrix(X)
  X_bar<-colWeightedMeans(X, w=prob, na.rm=T)
  X_bar<-X_bar[names(X_bar)%in%h]
  X_bar<-c(1,X_bar)
  names(X_bar)[1]<-h[1]
  return(X_bar)
}

Gini_hat<-function(X_bar,II,Q=c(0,0.25,0.5,0.75,1)){
  q=length(Q)
  Gini_hat<-matrix(rep(0,q),ncol=1,nrow=q)
  g<-apply(apply(II,2,'*',X_bar),2,sum)
  for (qq in 0:(q-2)){
    Gini_hat[qq+1]<-2*((Q[qq+2]-Q[qq+1])*Q[qq+1]+0.5*(Q[qq+2]-Q[qq+1])^2)
    Gini_hat[qq+1]<-Gini_hat[qq+1]-g[qq+1]
  }
  Gini_hat[q]=1-g[q]
  return(Gini_hat)
}

Ind_cov<-function(X0,X1,cov){
  temp<-X1
  temp[cov]=X0[cov]
  return(temp)
}

Ind_coef<-function(II0,II1,coef){
  temp<-II0
  temp[coef,]=II1[coef,]
  return(temp)
}

firstup <- function(x) {
  substr(x, 1, 1) <- toupper(substr(x, 1, 1))
  x
}

#Function to get time attribute
`%.%` <- function(o, a) attr(o, as.character(substitute(a)))

